module lambda0_optimizer_mod
USE functions_mod
USE parameters_mod
IMPLICIT NONE

CONTAINS
SUBROUTINE lambda0_optimizer(lambda0_out,gov_spend,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,&
                    types,partprod,rate,tr,population,avg_earnings,ss,survive,ss_cap_use,lambda1_use,lambdak_use)
INTEGER, intent(in)   :: sim_types(J,simulationsN),sim_shocks(J,simulationsN)
INTEGER, parameter ::  guessesN=30
DOUBLE PRECISION,    intent(in)   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),hbar,wage,survive(J)
DOUBLE PRECISION,    intent(in)   :: shocks(shocksN),types(typesN),sim_energy(J,simulationsN),ss_cap_use,lambda1_use,lambdak_use
DOUBLE PRECISION,    intent(in)   :: partprod,rate,tr,gov_spend,population(J),avg_earnings,ss(typesN,shocksN/2)
DOUBLE PRECISION,    intent(out)  :: lambda0_out
DOUBLE PRECISION :: diff_out,guess_l,guess_h,guess_first,test_l,test_h,test_1,test_2,guess_2,guesses(guessesN),test(guessesN)
INTEGER ::keep_going,counter,it,guessint_l,guessint_h

guess_l=.1D0
guess_h=1.1D0
guesses(1)=guess_l
test(1)=tax_differ_lambda0_finder(guesses(1),sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,&
    partprod,rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda1_use,lambdak_use)
guesses(guessesN)=guess_h
test(guessesN)=tax_differ_lambda0_finder(guesses(guessesN),sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,&
    partprod,rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda1_use,lambdak_use)

do it=2,guessesN-1
    guesses(it)=guesses(it-1)+(guess_h-guess_l)/(real(guessesN)-1.0D0)
    test(it)=tax_differ_lambda0_finder(guesses(it),sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,&
    partprod,rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda1_use,lambdak_use)
end do

if(test(guessesN)<0.0D0) then
    lambda0_out=guesses(guessesN)
elseif(test(1)>0.0D0) then
    lambda0_out=guesses(1)
else
    guessint_h=minloc(test,1,mask=test.ge.0)
    guess_h=guesses(guessint_h)
    guessint_l=maxloc(test,1,mask=test.le.0,back=.TRUE.)
    guess_l=guesses(guessint_l)
    guess_first=(guess_l+guess_h)/2.0D0

print*,guess_l,guess_h,guessint_l,guessint_H

diff_out=brent(guess_l, guess_first, guess_h,tax_differ_lambda0,.000000001D0,lambda0_out,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,&
                hbar,wage,shocks,types,partprod,rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda1_use,lambdak_use)

end if

end subroutine lambda0_optimizer

DOUBLE PRECISION FUNCTION tax_differ_lambda0_finder(lambda0_test,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,&
    partprod,rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda1_use,lambdak_use)
    INTEGER, intent(in)   :: sim_types(J,simulationsN),sim_shocks(J,simulationsN)
    DOUBLE PRECISION,    intent(in)   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),hbar,wage,survive(J),ss_cap_use
    DOUBLE PRECISION,    intent(in)   :: shocks(shocksN),types(typesN),sim_energy(J,simulationsN),lambda1_use,lambdak_use
    DOUBLE PRECISION,    intent(in)   :: partprod,rate,tr,lambda0_test,gov_spend,population(J),avg_earnings,ss(typesN,shocksN/2)
    DOUBLE PRECISION	::tax_rev,sim_taxespaid(J,simulationsN),labor_income,rate_use,labor_tax,capital_tax
    INTEGER ::it,s
    !! Taxes
    tax_rev=0.0D0
    sim_taxespaid=0.0D0
    do it=1,simulationsN
        do s=1,Jnr-1
            if(sim_capital(s,it)<0.0D0) then
                if(s==1) then
                    rate_use=rate/survive(1)
                else
                    rate_use=rate/survive(s-1)
                end if
            else
                rate_use=rate
            end if

            if (sim_labor(s,it)<hbar) then
                labor_income=wage*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod*sim_labor(s,it)
                sim_taxespaid(s,it)=all_taxer(rate_use*(sim_capital(s,it)+tr),&
                    labor_income,1,&
                    ss_cap_use,lambda0_test,lambda1_use,lambdak_use)+taue*sim_energy(s,it)
            else
                labor_income=wage*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*sim_labor(s,it)
                sim_taxespaid(s,it)=all_taxer(rate_use*(sim_capital(s,it)+tr),&
                    labor_income,1,ss_cap_use,lambda0_test,lambda1_use,lambdak_use)&
                    +taue*sim_energy(s,it)
            end if
            tax_rev=tax_rev+sim_taxespaid(s,it)*population(s)/simulationsN
        end do
        do s=Jnr,J
                sim_taxespaid(s,it)=all_taxer(rate*(sim_capital(s,it)+tr),0.0D0,0,ss_cap_use,lambda0_test,lambda1_use,lambdak_use)+&
                    taue*sim_energy(s,it)


            tax_rev=tax_rev+sim_taxespaid(s,it)*population(s)/simulationsN
        end do
    end do
    tax_differ_lambda0_finder=(gov_spend-tax_rev)
end FUNCTION tax_differ_lambda0_finder


!This is objective function to find lambda0
DOUBLE PRECISION FUNCTION tax_differ_lambda0(lambda0_test,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,&
    partprod,rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda1_use,lambdak_use)
    INTEGER, intent(in)   :: sim_types(J,simulationsN),sim_shocks(J,simulationsN)
    DOUBLE PRECISION,    intent(in)   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),hbar,wage,survive(J),ss_cap_use
    DOUBLE PRECISION,    intent(in)   :: shocks(shocksN),types(typesN),sim_energy(J,simulationsN),lambda1_use,lambdak_use
    DOUBLE PRECISION,    intent(in)   :: partprod,rate,tr,lambda0_test,gov_spend,population(J),avg_earnings,ss(typesN,shocksN/2)
    DOUBLE PRECISION	::tax_rev,sim_taxespaid(J,simulationsN),labor_income,rate_use,labor_tax,capital_tax
    INTEGER ::it,s
    !! Taxes
    tax_rev=0.0D0
    sim_taxespaid=0.0D0
    do it=1,simulationsN
        do s=1,Jnr-1
            if(sim_capital(s,it)<0.0D0) then
                if(s==1) then
                    rate_use=rate/survive(1)
                else
                    rate_use=rate/survive(s-1)
                end if
            else
                rate_use=rate
            end if

            if (sim_labor(s,it)<hbar) then
                labor_income=wage*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod*sim_labor(s,it)
                sim_taxespaid(s,it)=all_taxer(rate_use*(sim_capital(s,it)+tr),&
                    labor_income,1,&
                    ss_cap_use,lambda0_test,lambda1_use,lambdak_use)+taue*sim_energy(s,it)
            else
                labor_income=wage*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*sim_labor(s,it)
                sim_taxespaid(s,it)=all_taxer(rate_use*(sim_capital(s,it)+tr),&
                    labor_income,1,ss_cap_use,lambda0_test,lambda1_use,lambdak_use)&
                    +taue*sim_energy(s,it)
            end if
            tax_rev=tax_rev+sim_taxespaid(s,it)*population(s)/simulationsN
        end do
        do s=Jnr,J
                sim_taxespaid(s,it)=all_taxer(rate*(sim_capital(s,it)+tr),0.0D0,0,ss_cap_use,lambda0_test,lambda1_use,lambdak_use)+&
                    taue*sim_energy(s,it)


            tax_rev=tax_rev+sim_taxespaid(s,it)*population(s)/simulationsN
        end do
    end do
    tax_differ_lambda0=abs(gov_spend-tax_rev)
end FUNCTION tax_differ_lambda0



!This is objective function to find lambda0
DOUBLE PRECISION FUNCTION tax_differ_lambda01(lambda0_test,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,&
    partprod,rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda1_use,lambdak_use)
    INTEGER, intent(in)   :: sim_types(J,simulationsN),sim_shocks(J,simulationsN)
    DOUBLE PRECISION,    intent(in)   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),hbar,wage,survive(J),ss_cap_use
    DOUBLE PRECISION,    intent(in)   :: shocks(shocksN),types(typesN),sim_energy(J,simulationsN),lambda1_use,lambdak_use
    DOUBLE PRECISION,    intent(in)   :: partprod,rate,tr,lambda0_test,gov_spend,population(J),avg_earnings,ss(typesN,shocksN/2)
    DOUBLE PRECISION	::tax_rev,sim_taxespaid(J,simulationsN),labor_income,rate_use,labor_tax,capital_tax
    INTEGER ::it,s
    !! Taxes
    tax_rev=0.0D0
    sim_taxespaid=0.0D0
    do it=1,simulationsN
        do s=1,Jnr-1
            if(sim_capital(s,it)<0.0D0) then
                if(s==1) then
                    rate_use=rate/survive(1)
                else
                    rate_use=rate/survive(s-1)
                end if
            else
                rate_use=rate
            end if

            if (sim_labor(s,it)<hbar) then
                labor_income=wage*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod*sim_labor(s,it)
                sim_taxespaid(s,it)=all_taxer(rate_use*(sim_capital(s,it)+tr),&
                    labor_income,1,&
                    ss_cap_use,lambda0_test,lambda1_use,lambdak_use)+taue*sim_energy(s,it)
            else
                labor_income=wage*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*sim_labor(s,it)
                sim_taxespaid(s,it)=all_taxer(rate_use*(sim_capital(s,it)+tr),&
                    labor_income,1,ss_cap_use,lambda0_test,lambda1_use,lambdak_use)&
                    +taue*sim_energy(s,it)
            end if
            tax_rev=tax_rev+sim_taxespaid(s,it)*population(s)/simulationsN
        end do
        do s=Jnr,J
                sim_taxespaid(s,it)=all_taxer(rate*(sim_capital(s,it)+tr),0.0D0,0,ss_cap_use,lambda0_test,lambda1_use,lambdak_use)+&
                    taue*sim_energy(s,it)


            tax_rev=tax_rev+sim_taxespaid(s,it)*population(s)/simulationsN
        end do
    end do
    tax_differ_lambda01=abs(gov_spend-tax_rev)
    print*,"lambda out",lambda0_test,gov_spend,tax_rev,tax_differ_lambda01
end FUNCTION tax_differ_lambda01







SUBROUTINE lambdak_optimizer(lambdak_out,gov_spend,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,&
                    types,partprod,rate,tr,population,avg_earnings,ss,survive,ss_cap_use,lambda0_use,lambda1_use,lambdak_use)
INTEGER, intent(in)   :: sim_types(J,simulationsN),sim_shocks(J,simulationsN)
DOUBLE PRECISION,    intent(in)   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),hbar,wage,survive(J),lambdak_use
DOUBLE PRECISION,    intent(in)   :: shocks(shocksN),types(typesN),sim_energy(J,simulationsN),ss_cap_use,lambda0_use,lambda1_use
DOUBLE PRECISION,    intent(in)   :: partprod,rate,tr,gov_spend,population(J),avg_earnings,ss(typesN,shocksN/2)
DOUBLE PRECISION,    intent(out)  :: lambdak_out
DOUBLE PRECISION :: diff_out



diff_out=brent_tax(-0.30D0, lambdak_use, 0.80D0,tax_differ_lambdak,.000000001D0,lambdak_out,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,&
                hbar,wage,shocks,types,partprod,rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda0_use,lambda1_use)
end subroutine lambdak_optimizer


!This is objective function to find lambda0
DOUBLE PRECISION FUNCTION tax_differ_lambdak(lambdak_test,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,&
    partprod,rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda0_use,lambda1_use)
    INTEGER, intent(in)   :: sim_types(J,simulationsN),sim_shocks(J,simulationsN)
    DOUBLE PRECISION,    intent(in)   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),hbar,wage,survive(J)
    DOUBLE PRECISION,    intent(in)   :: shocks(shocksN),types(typesN),sim_energy(J,simulationsN),ss_cap_use,lambda0_use,lambda1_use
    DOUBLE PRECISION,    intent(in)   :: partprod,rate,tr,lambdak_test,gov_spend,population(J),avg_earnings,ss
    DOUBLE PRECISION	::tax_rev,sim_taxespaid(J,simulationsN),labor_income,rate_use,labor_tax,capital_tax
    INTEGER ::it,s
    !! Taxes
    tax_rev=0.0D0
    sim_taxespaid=0.0D0
    do it=1,simulationsN
        do s=1,Jnr-1
            if(sim_capital(s,it)<0.0D0) then
                if(s==1) then
                    rate_use=rate/survive(1)
                else
                    rate_use=rate/survive(s-1)
                end if
            else
                rate_use=rate
            end if

            if (sim_labor(s,it)<hbar) then
                labor_income=wage*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod*sim_labor(s,it)
                sim_taxespaid(s,it)=all_taxer_k(rate_use*(sim_capital(s,it)+tr),&
                    labor_income,1,&
                    lambdak_test,ss_cap_use,lambda0_use,lambda1_use)+taue*sim_energy(s,it)
            else
                labor_income=wage*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*sim_labor(s,it)
                sim_taxespaid(s,it)=all_taxer_k(rate_use*(sim_capital(s,it)+tr),&
                    labor_income,1,lambdak_test,ss_cap_use,lambda0_use,lambda1_use)&
                    +taue*sim_energy(s,it)
            end if
            tax_rev=tax_rev+sim_taxespaid(s,it)*population(s)/simulationsN
        end do
        do s=Jnr,J
            sim_taxespaid(s,it)=all_taxer_k(rate*(sim_capital(s,it)+tr),0.0D0,0,lambdak_test,ss_cap_use,lambda0_use,lambda1_use)+&
            taue*sim_energy(s,it)

            tax_rev=tax_rev+sim_taxespaid(s,it)*population(s)/simulationsN
        end do
    end do
    tax_differ_lambdak=abs(gov_spend-tax_rev)
end FUNCTION tax_differ_lambdak










SUBROUTINE rebate_optimizer(rebate_out,gov_spend,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,&
                    types,partprod,rate,tr,population,avg_earnings,ss,survive,ss_cap_use,lambda0_use,lambda1_use,lambdak_use)
INTEGER, intent(in)   :: sim_types(J,simulationsN),sim_shocks(J,simulationsN)
DOUBLE PRECISION,    intent(in)   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),hbar,wage,survive(J)
DOUBLE PRECISION,    intent(in)   :: shocks(shocksN),types(typesN),sim_energy(J,simulationsN),ss_cap_use,lambda0_use,lambda1_use,lambdak_use
DOUBLE PRECISION,    intent(in)   :: partprod,rate,tr,gov_spend,population(J),avg_earnings,ss(typesN,shocksN/2)
DOUBLE PRECISION,    intent(out)  :: rebate_out
DOUBLE PRECISION :: diff_out



diff_out=brent_rebate(0.0D0, rebate_level(year_calculating), 0.50D0,tax_differ_rebate,.000000001D0,rebate_out,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,&
                hbar,wage,shocks,types,partprod,rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda0_use,lambda1_use,lambdak_use)
end subroutine rebate_optimizer


!This is objective function to find lambda0
DOUBLE PRECISION FUNCTION tax_differ_rebate(rebate_test,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,&
    partprod,rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda0_use,lambda1_use,lambdak_use)
    INTEGER, intent(in)   :: sim_types(J,simulationsN),sim_shocks(J,simulationsN)
    DOUBLE PRECISION,    intent(in)   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),hbar,wage,survive(J)
    DOUBLE PRECISION,    intent(in)   :: shocks(shocksN),types(typesN),sim_energy(J,simulationsN),ss_cap_use,lambda0_use,lambda1_use,lambdak_use
    DOUBLE PRECISION,    intent(in)   :: partprod,rate,tr,rebate_test,gov_spend,population(J),avg_earnings,ss(typesN,shocksN/2)
    DOUBLE PRECISION	::tax_rev,sim_taxespaid(J,simulationsN),labor_income,rate_use,labor_tax,capital_tax
    INTEGER ::it,s
    !! Taxes
    tax_rev=0.0D0
    sim_taxespaid=0.0D0
    do it=1,simulationsN
        do s=1,Jnr-1
            if(sim_capital(s,it)<0.0D0) then
                if(s==1) then
                    rate_use=rate/survive(1)
                else
                    rate_use=rate/survive(s-1)
                end if
            else
                rate_use=rate
            end if

            if (sim_labor(s,it)<hbar) then
                labor_income=wage*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod*sim_labor(s,it)
                sim_taxespaid(s,it)=all_taxer_rebate(rate_use*(sim_capital(s,it)+tr),&
                    labor_income,1,rebate_test,ss_cap_use,lambda0_use,lambda1_use,lambdak_use)+taue*sim_energy(s,it)
            else
                labor_income=wage*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*sim_labor(s,it)
                sim_taxespaid(s,it)=all_taxer_rebate(rate_use*(sim_capital(s,it)+tr),&
                    labor_income,1,rebate_test,ss_cap_use,lambda0_use,lambda1_use,lambdak_use)&
                    +taue*sim_energy(s,it)
            end if
            tax_rev=tax_rev+sim_taxespaid(s,it)*population(s)/simulationsN
        end do
        do s=Jnr,J
            sim_taxespaid(s,it)=all_taxer_rebate(rate*(sim_capital(s,it)+tr),0.0D0,0,rebate_test,ss_cap_use,lambda0_use,lambda1_use,lambdak_use)+&
                taue*sim_energy(s,it)
            tax_rev=tax_rev+sim_taxespaid(s,it)*population(s)/simulationsN
        end do
    end do
    tax_differ_rebate=abs(gov_spend-tax_rev)
end FUNCTION tax_differ_rebate




!This is function to find out difference
DOUBLE PRECISION FUNCTION tax_differ2(lambda0_test,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,&
    partprod,rate,tr,gov_spend,population,avg_earnings,ss,lambdak_test,survive,rebate_test,ss_cap_use,lambda1_use)
    INTEGER, intent(in)   :: sim_types(J,simulationsN),sim_shocks(J,simulationsN)
    DOUBLE PRECISION,    intent(in)   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),hbar,wage,survive(J),rebate_test
    DOUBLE PRECISION,    intent(in)   :: shocks(shocksN),types(typesN),sim_energy(J,simulationsN),lambdak_test,ss_cap_use,lambda1_use
    DOUBLE PRECISION,    intent(in)   :: partprod,rate,tr,lambda0_test,gov_spend,population(J),avg_earnings,ss(typesN,shocksN/2)
    DOUBLE PRECISION	::tax_rev,sim_taxespaid(J,simulationsN),labor_income,rate_use,labor_tax,capital_tax
    INTEGER ::it,s
    !! Taxes
    tax_rev=0.0D0
    sim_taxespaid=0.0D0
    do it=1,simulationsN
        do s=1,Jnr-1
            if(sim_capital(s,it)<0.0D0) then
                if(s==1) then
                    rate_use=rate/survive(1)
                else
                    rate_use=rate/survive(s-1)
                end if
            else
                rate_use=rate
            end if

            if (sim_labor(s,it)<hbar) then
                labor_income=wage*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod*sim_labor(s,it)
                sim_taxespaid(s,it)=all_taxer_b(rate_use*(sim_capital(s,it)+tr),&
                    labor_income,1,lambda0_test,lambdak_test,rebate_test,ss_cap_use,lambda1_use)+taue*sim_energy(s,it)
            else
                labor_income=wage*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*sim_labor(s,it)
                sim_taxespaid(s,it)=all_taxer_b(rate_use*(sim_capital(s,it)+tr),&
                    labor_income,1,lambda0_test,lambdak_test,rebate_test,ss_cap_use,lambda1_use)&
                    +taue*sim_energy(s,it)
            end if
            tax_rev=tax_rev+sim_taxespaid(s,it)*population(s)/simulationsN
        end do
        do s=Jnr,J
                sim_taxespaid(s,it)=all_taxer_b(rate*(sim_capital(s,it)+tr),0.0D0,0,lambda0_test,lambdak_test,rebate_test,&
                    ss_cap_use,lambda1_use)+&
                    taue*sim_energy(s,it)


            tax_rev=tax_rev+sim_taxespaid(s,it)*population(s)/simulationsN
        end do
    end do
    tax_differ2=(gov_spend-tax_rev)
end FUNCTION tax_differ2




      DOUBLE PRECISION FUNCTION brent(ax,bx,cx,func,tol,xmin,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,&
                    shocks,types,partprod,rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda1_use,lambdak_use)
USE nrtype; USE nrutil, ONLY : nrerror
IMPLICIT NONE
    INTEGER, intent(in)   :: sim_types(J,simulationsN),sim_shocks(J,simulationsN)
    DOUBLE PRECISION, INTENT(IN)    :: ax,bx,cx,tol
    DOUBLE PRECISION,    intent(in)   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),hbar,wage,survive(J)
    DOUBLE PRECISION,    intent(in)   :: shocks(shocksN),types(typesN),sim_energy(J,simulationsN),ss_cap_use,lambda1_use,lambdak_use
    DOUBLE PRECISION,    intent(in)   :: partprod,rate,tr,gov_spend,population(J),avg_earnings,ss(typesN,shocksN/2)
    DOUBLE PRECISION, INTENT(OUT)   :: xmin
INTEGER, PARAMETER          :: ITMAX=500
DOUBLE PRECISION, PARAMETER ::CGOLD=0.3819660D0,ZEPS=1.0e-3_sp*epsilon(ax)
INTEGER                     ::iter
DOUBLE PRECISION            :: a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
DOUBLE PRECISION :: func
a=min(ax,cx)
b=max(ax,cx)
v=bx
w=v
x=v
e=0.0
fx=func(x,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,partprod,&
                    rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda1_use,lambdak_use)
fv=fx
fw=fx
do iter=1,ITMAX
    xm=0.5D0*(a+b)
    tol1=tol*abs(x)+ZEPS
    tol2=2.0D0*tol1
    if (abs(x-xm) <= (tol2-0.5D0*(b-a))) then
        xmin=x
        brent=fx
        RETURN
    end if
    if (abs(e) > tol1) then
        r=(x-w)*(fx-fv)
        q=(x-v)*(fx-fw)
        p=(x-v)*q-(x-w)*r
        q=2.0D0*(q-r)
        if (q > 0.0D0) p=-p
        q=abs(q)
        etemp=e
        e=d
        if (abs(p) >= abs(0.5D0*q*etemp) .or. &
            p <= q*(a-x) .or. p >= q*(b-x)) then
            e=merge(a-x,b-x, x >= xm )
            d=CGOLD*e
        else
            d=p/q
            u=x+d
            if (u-a < tol2 .or. b-u < tol2) d=sign(tol1,xm-x)
        end if
    else
        e=merge(a-x,b-x, x >= xm )
        d=CGOLD*e
    end if
    u=merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
    fu=func(u,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,partprod,&
                    rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda1_use,lambdak_use)
    if (fu <= fx) then
        if (u >= x) then
            a=x
        else
            b=x
        end if
        call shft(v,w,x,u)
        call shft(fv,fw,fx,fu)
    else
        if (u < x) then
            a=u
        else
            b=u
        end if
        if (fu <= fw .or. w == x) then
            v=w
            fv=fw
            w=u
            fw=fu
        else if (fu <= fv .or. v == x .or. v == w) then
            v=u
            fv=fu
        end if
    end if
end do
call nrerror('brent lambda: exceed maximum iterations')
CONTAINS
!BL
SUBROUTINE shft(a,b,c,d)
    DOUBLE PRECISION, INTENT(OUT)   :: a
    DOUBLE PRECISION, INTENT(INOUT) :: b,c
    DOUBLE PRECISION, INTENT(IN)    :: d
a=b
b=c
c=d
END SUBROUTINE shft
END FUNCTION brent





      DOUBLE PRECISION FUNCTION brent_rebate(ax,bx,cx,func,tol,xmin,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,&
                    shocks,types,partprod,rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda0_use,lambda1_use,lambdak_use)
USE nrtype; USE nrutil, ONLY : nrerror
IMPLICIT NONE
    INTEGER, intent(in)   :: sim_types(J,simulationsN),sim_shocks(J,simulationsN)
    DOUBLE PRECISION, INTENT(IN)    :: ax,bx,cx,tol
    DOUBLE PRECISION,    intent(in)   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),hbar,wage,survive(J)
    DOUBLE PRECISION,    intent(in)   :: shocks(shocksN),types(typesN),sim_energy(J,simulationsN),ss_cap_use,lambda1_use,lambdak_use
    DOUBLE PRECISION,    intent(in)   :: partprod,rate,tr,gov_spend,population(J),avg_earnings,ss(typesN,shocksN/2),lambda0_use
    DOUBLE PRECISION, INTENT(OUT)   :: xmin
INTEGER, PARAMETER          :: ITMAX=500
DOUBLE PRECISION, PARAMETER ::CGOLD=0.3819660D0,ZEPS=1.0e-3_sp*epsilon(ax)
INTEGER                     ::iter
DOUBLE PRECISION            :: a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
DOUBLE PRECISION :: func
a=min(ax,cx)
b=max(ax,cx)
v=bx
w=v
x=v
e=0.0
fx=func(x,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,partprod,&
                    rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda0_use,lambda1_use,lambdak_use)
fv=fx
fw=fx
do iter=1,ITMAX
    xm=0.5D0*(a+b)
    tol1=tol*abs(x)+ZEPS
    tol2=2.0D0*tol1
    if (abs(x-xm) <= (tol2-0.5D0*(b-a))) then
        xmin=x
        brent_rebate=fx
        RETURN
    end if
    if (abs(e) > tol1) then
        r=(x-w)*(fx-fv)
        q=(x-v)*(fx-fw)
        p=(x-v)*q-(x-w)*r
        q=2.0D0*(q-r)
        if (q > 0.0D0) p=-p
        q=abs(q)
        etemp=e
        e=d
        if (abs(p) >= abs(0.5D0*q*etemp) .or. &
            p <= q*(a-x) .or. p >= q*(b-x)) then
            e=merge(a-x,b-x, x >= xm )
            d=CGOLD*e
        else
            d=p/q
            u=x+d
            if (u-a < tol2 .or. b-u < tol2) d=sign(tol1,xm-x)
        end if
    else
        e=merge(a-x,b-x, x >= xm )
        d=CGOLD*e
    end if
    u=merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
    fu=func(u,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,partprod,&
                    rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda0_use,lambda1_use,lambdak_use)
    if (fu <= fx) then
        if (u >= x) then
            a=x
        else
            b=x
        end if
        call shft(v,w,x,u)
        call shft(fv,fw,fx,fu)
    else
        if (u < x) then
            a=u
        else
            b=u
        end if
        if (fu <= fw .or. w == x) then
            v=w
            fv=fw
            w=u
            fw=fu
        else if (fu <= fv .or. v == x .or. v == w) then
            v=u
            fv=fu
        end if
    end if
end do
call nrerror('brent_rebate lambda: exceed maximum iterations')
CONTAINS
!BL
SUBROUTINE shft(a,b,c,d)
    DOUBLE PRECISION, INTENT(OUT)   :: a
    DOUBLE PRECISION, INTENT(INOUT) :: b,c
    DOUBLE PRECISION, INTENT(IN)    :: d
a=b
b=c
c=d
END SUBROUTINE shft
END FUNCTION brent_rebate



      DOUBLE PRECISION FUNCTION brent_tax(ax,bx,cx,func,tol,xmin,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,&
                    shocks,types,partprod,rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda0_use,lambda1_use)
USE nrtype; USE nrutil, ONLY : nrerror
IMPLICIT NONE
    INTEGER, intent(in)   :: sim_types(J,simulationsN),sim_shocks(J,simulationsN)
    DOUBLE PRECISION, INTENT(IN)    :: ax,bx,cx,tol
    DOUBLE PRECISION,    intent(in)   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),hbar,wage,survive(J)
    DOUBLE PRECISION,    intent(in)   :: shocks(shocksN),types(typesN),sim_energy(J,simulationsN),ss_cap_use,lambda0_use,lambda1_use
    DOUBLE PRECISION,    intent(in)   :: partprod,rate,tr,gov_spend,population(J),avg_earnings,ss(typesN,shocksN/2)
    DOUBLE PRECISION, INTENT(OUT)   :: xmin
INTEGER, PARAMETER          :: ITMAX=500
DOUBLE PRECISION, PARAMETER ::CGOLD=0.3819660D0,ZEPS=1.0e-3_sp*epsilon(ax)
INTEGER                     ::iter
DOUBLE PRECISION            :: a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
DOUBLE PRECISION :: func
a=min(ax,cx)
b=max(ax,cx)
v=bx
w=v
x=v
e=0.0
fx=func(x,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,partprod,&
                    rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda0_use,lambda1_use)
fv=fx
fw=fx
do iter=1,ITMAX
    xm=0.5D0*(a+b)
    tol1=tol*abs(x)+ZEPS
    tol2=2.0D0*tol1
    if (abs(x-xm) <= (tol2-0.5D0*(b-a))) then
        xmin=x
        brent_tax=fx
        RETURN
    end if
    if (abs(e) > tol1) then
        r=(x-w)*(fx-fv)
        q=(x-v)*(fx-fw)
        p=(x-v)*q-(x-w)*r
        q=2.0D0*(q-r)
        if (q > 0.0D0) p=-p
        q=abs(q)
        etemp=e
        e=d
        if (abs(p) >= abs(0.5D0*q*etemp) .or. &
            p <= q*(a-x) .or. p >= q*(b-x)) then
            e=merge(a-x,b-x, x >= xm )
            d=CGOLD*e
        else
            d=p/q
            u=x+d
            if (u-a < tol2 .or. b-u < tol2) d=sign(tol1,xm-x)
        end if
    else
        e=merge(a-x,b-x, x >= xm )
        d=CGOLD*e
    end if
    u=merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
    fu=func(u,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,wage,shocks,types,partprod,&
                    rate,tr,gov_spend,population,avg_earnings,ss,survive,ss_cap_use,lambda0_use,lambda1_use)
    if (fu <= fx) then
        if (u >= x) then
            a=x
        else
            b=x
        end if
        call shft(v,w,x,u)
        call shft(fv,fw,fx,fu)
    else
        if (u < x) then
            a=u
        else
            b=u
        end if
        if (fu <= fw .or. w == x) then
            v=w
            fv=fw
            w=u
            fw=fu
        else if (fu <= fv .or. v == x .or. v == w) then
            v=u
            fv=fu
        end if
    end if
end do
call nrerror('brent_tax lambda: exceed maximum iterations')
CONTAINS
!BL
SUBROUTINE shft(a,b,c,d)
    DOUBLE PRECISION, INTENT(OUT)   :: a
    DOUBLE PRECISION, INTENT(INOUT) :: b,c
    DOUBLE PRECISION, INTENT(IN)    :: d
a=b
b=c
c=d
END SUBROUTINE shft
END FUNCTION brent_tax

end module lambda0_optimizer_mod
